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Abstract 

We investigate the dynamics of a binary mixture Lennard-Jones system of different system sizes 
with respect to the importance of the properties of the underlying potential energy landscape 
(PEL). We show that the dynamics of small systems can be very well described within the con- 
tinuous time random walk formalism, which is determined solely by PEL parameters. Finite size 
analysis shows that the diffusivity of large and small systems are very similar. This suggests that 
the PEL parameters of the small system also determine the local dynamics in large systems. The 
structural relaxation time, however, displays significant finite size effects. Furthermore, using a 
non-equilibrium configuration of a large system, we find that causal connections exist between 
close-by regions of the system. These findings can be described by the coupled landscape model for 
which a macroscopic system is described by a superposition of elementary systems each described 
by its PEL. A minimum coupling is introduced which accounts for the finite size behavior. The 
coupling strength, as the single adjustable parameter, becomes smaller closer to the glass transition. 



I. INTRODUCTION 



Understanding the complex dynamics of glassy materials such as supercooled liquids 
is still a highly controversial problem. In contrast to regular liquids the single particle 
dynamics in the supercooled state become spatially and temporally correlated, which leads to 
pronounced dynamical heterogeneities. Dynamical heterogeneities are a universal property 
of glass-forming systems and are responsible for several phenomena like non-exponential 
decay of response functions and the violation of the Stokes-Einstein relation. 

In small model systems Vogel and coworkers pQ could show that particle rearrangements 
are typically localized and that their number does not depend on temperature. These find- 
ings could also be validated for large systems by Keys et al. |2j. The authors additionally 
related the propagation of motion to dynamical subunits which occur in a facilitation-like 
manner, meaning that excitations on one scale facilitate dynamics of neighboring excitations 
thereby creating excitations on larger scales. The characterization of facilitated dynamics 
in molecular dynamics (MD) simulations was attempted by different authors. For exam- 
ple Candelier and coworkers (3j H] proved the aggregation of cage-breaking processes into 
avalanches in a granular system. Presently discussed models range between the old idea 
of Adam and Gibbs j5] of a qualitative decomposition of a sample into elementary subsys- 
tems and the facilitation approach, where all complexity emerges by more or less complex 
coupling rules. For example the kinetically constrained spin models were studied exten- 
sively by different authors [El [7]. Dynamical facilitation can alternatively be formulated 
as a coupling process. Rehwald et al. [8] discovered from studying finite size effects, that 
these coupling processes determine structural relaxation properties. In shear experiments it 
was demonstrated that the range of coupling interactions, manifested by correlated particle 
displacements due to stress, covers just some particle diameters [HI EI]. Other groups report 
of long ranged elastic interactions in glycerol measured by dielectric relaxation [TT] . 

A different approach, relating the dramatic slowing down of the dynamics at the glass 
transition to the potential energy landscape (PEL), was introduced by Goldstein [12J. He 
stated that at sufficiently low temperatures the system resides near local minima, so called 
inherent structures (IS), of the potential. While many computer simulations confirmed that 
the thermodynamics are indeed determined by the distribution of ISs, the connection with 
the dynamics could be revealed much later. For small systems it was possible to predict 
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the dynamics from parameters of the underlying potential energy landscape (PEL) [T3TfT6] 
with the concept of metabasins. The dynamics between metabasins are closely related to 
the trap model [T7J [TS]. Many features of the complex glassy dynamics could be derived 
from PEL properties. 

The situation changes when studying large systems: the PEL is no longer suitable for 
the prediction of the dynamics. Here we show that the PEL formalism can be extended 
to large systems by decomposing the system into elementary subsystems each of which is 
described by its own PEL. In contrast to kinetically constrained spin models the elementary 
subunits in this approach are complex systems, already reflecting macroscopic properties for 
thermodynamic observables pj5] or self diffusion properties. Since thermodynamic proper- 
ties do not change significantly when enlarging the system above a certain minimum size 
[19], and rearrangements are localized within a temperature independent size, we can be 
sure to capture the important relaxation properties in a single elementary subsystem. Ad- 
ditional fluctuations and interactions can then be included into the coupling rules between 
the subsystems. In terms of the complex nature of the elementary subsystem our approach 
resembles the mosaic approach [20] where the sample is decomposed into a mosaic of aperi- 
odic crystals, so called entropic droplets. Nevertheless, in the mosaic approach all coupling 
processes are captured by the distribution of free energy barrier heights. 

This manuscript is organized as follows: First we show how the minimum system can 
be described by the PEL and the continuous time random walk approach (CTRW), both 
analytically and numerically. After discussing changes of the CTRW when approaching 
larger system sizes we give a short summary of the physical scenario and motivate the 
presence of coupling effects. After discussing some technical questions concerning the used 
variables for the comparison between model and MD, we explicitly identify the presence 
of coupling processes via appropriate simulations. These findings motivate the coupled 
landscape model (CLM) which will be presented in the last section. It is shown that to a 
very good approximation the finite size effects of dynamic observables are fully captured 
by the CLM. Finally we figure out in how far the CLM is connected to recently discussed 
models. 
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TABLE I. Thermodynamic and dynamic PEL parameters for the binary mixture Lennard-Jones 
system with N = 65 and eo = 0. 

II. THE MINIMUM SYSTEM 

Throughout this paper we compare the model results with the binary mixture Kob- 
Andersen Lennard-Jones (BMLJ) system [21] applying periodic boundary conditions. Due 
to the small system size we have used a slightly shorter cutoff of r c = 1.8 |14j . All MD 
simulations have been performed in the NVT ensemble using a Nose-Hoover thermostat. 

The PEL of a binary mixture of Lennard-Jones particles applying different system sizes 
was studied extensively [22]. One key result of Doliwa and Heuer [TH] is the connection 
between states of the PEL, so called metabasins (MB), and the dynamics of the system. 
An MB can be constructed from a given trajectory of inherent structures by removing the 
complete forward-backward correlations between them. If the system is small enough, the 
mean waiting time (r(e)), during which the system resides in a MB, strongly depends on the 
energy e of the specific MB. Of course, the system cannot be too small because otherwise 
massive finite size effects set in regarding the thermodynamics [Hj. In the same work it 
has been shown, that a system size of iV m in = 65 (BMLJ65) particles is a good choice for a 
minimum system. This minimum length scale of approx. 65 particles does not display any 
significant temperature dependence. 

In the next step one has to determine how the specific properties of the potential energy 
landscape determine the dynamic properties of the liquid. Many properties will resemble 
the simple trapmodel [T7J [23]. Here we summarize the most important results. The cor- 
responding parameters, determined from appropriate simulations (22] are listed in Tab. |TJ 
Due to improved simulation data they slightly differ from those given in |22j . 

The shape of the density of states (DOS) G(e) turned out to be close to a Gaussian 
G(e) ~ exp[— (e — eo)/2a 2 ] with an additional cutoff below e cut . For low e the G(e) decays 
slightly faster than expected from a Gaussian. For reasons of simplicity this is modeled by a 
cutoff energy e cut and an additional factor exp[— (e — e cut Y]. For analytical results this factor 
will be neglected. In the trap model the escape out of a MB of energy e can be described by 
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a simple rate process with rate T(e), given by T(e) = Foexp[/3e]. For the BMLJ65 system a 
slightly more complex energy dependence is observed [22] : 



c 



r(e) _ RV eX P [M/ 3 + «^entro)(e - e c )] , e < 

Fo (l ,e>e c 

where T defines the the overall time scale. The escape from a MB of energy e is a multi-step 
process. The energy of the first IS after having completed the escape process is denoted by 
e c . Due to percolation arguments e c is independent of e [22J. The energy at the final barrier 
is given by Vo + e c , i.e. Vo denotes the height of the last barrier crossed. The relaxation 
is solid-like (activated) for e < e c , whereas it is liquid-like otherwise [23]. For e < e c an 
additional entropic term, involving the factor nk en t IO , has to be taken into account, where 
^entro = (^o — e c )/cx 2 . It reflects that the number of escape paths from a MB increases 
exponentially with decreasing MB energy. The limit k = 1 can be rationalized in a simple 
PEL model |25j. For reasons of simplicity we chose the energy scale such that eo = from 
now on. In the general case k < 1 plays the role of an empirical factor of the order of one. 
If the investigated system contains M independent subsystems one gets A = 1/M. Some 
additional implications of the choice A < 1 will be discussed in the appendix. 

As reported in J2B], the width of the waiting time distribution at fixed energy e, expressed 
via 5(e) = (T 2 (e)) / (j(e)) 2 — 1, show deviations from an exponential distribution which 
disagrees with a pure rate process. This is, however, expected if the total system is a 
superposition of subsystems [22]. This broadening results from the fact that the total energy 
e can be decomposed in different ways. Formally, this broadening can be expressed by 
a distribution ^(7, e) of jump rates 7 at fixed energy e. Here we approximate the 7- 
dependence of 07(7, e) by a log-normal distribution with variance <r 7 and mean value \i{e). 
Instead of an exponential distribution one gets u>(r,e) = f djzu^'j, e)7exp[— jt] as waiting 
time distribution for states of energy e. The moments of this distribution are related to the 
moments of 7 via 

In the MD simulation only (T n (e)) can be calculated. Therefore one can estimate 07 by the 
relation 



a 2 = ln 



2 (r(e))l 



(3) 
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FIG. 1. Comparison of the absolute values of D and r Q from the MB trajectory with corresponding 
expressions from the CTRW PEL description. Included is the diagonal. 

We calculate the right hand side of equation ^ from MB trajectories and find that o~ 7 
is approx. 1.0, corresponding to a broadening of one order of magnitude. It is slightly 
temperature dependent and can to a good approximation be assumed to be constant in the 
relevant energy range. The distribution 07(7, e) is very narrow compared to the equilibrium 
distribution of rates p(logT) (standard deviation at T = 0.5 is approx. 4.4). We note in 
passing, that for a-p > the energy of a state is no longer solely sufficient to describe the 
specific rate 7 of a MB. 

As shown in [27J the dynamics can be described as a CTRW in configuration space. 
Therefore the most relevant transport coefficients like D or r a can be calculated analytically 
by solving the Gaussian integrals J T n (e)p eq (e)de [28]. Beyond the high temperature limit 
(no influence of the crossover energy) results for A = k = 1 can be found in literature like 
the well-known quadratic behavior r a oc exp[a 2 (/3 — A)) 2 ] with (3 = —k cntm [2J [22j 122]. For 
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general A < 1 the T-dependence of D and r a is not purely quadratic (see Appendix |A|) . 
However, for the Stokes-Einstein relation one again obtains a quadratic dependence 

In [Dr a ] ~ In [(i^ 1 ) (r)] =A V (/3 + ^ entro ) 2 . (4) 

In the low temperature limit one gets Dr a ~ r| with £ = 2A/(2 + A) which for the present 
case gives £ = 0.4. 

In what follows we compare the CTRW/PEL predictions with the actual MD data in 
the temperature interval [0.45,0.6]. The lowest temperature is close to the mode-coupling 
temperature T c [TJ. More specifically we compare three different types of trajectories: (i) 
continuous trajectory from standard MD simulation, (ii) hopping trajectory between the 
MB configurations resulting form the same MD simulation, (iii) CTRW trajectory, based on 
the waiting times which are generated from the above PEL approach. 

In a first step we determine the diffusion constant D. By construction (i) and (ii) will 
result in identical values of D. For the CTRW trajectory one directly obtains D = a 2 /6N (r) 
where a 2 is a weakly temperature dependent value and can be determined independently 
[50] . As shown in Fig. [I] the temperature dependent first moment of the waiting time 
distribution perfectly reflects the T-dependence of D(T). Of course, strictly speaking this 
comparison is just a consistency check of the CTRW/PEL approach. 

In the second step we compare the structural relaxation times r a . In the CTRW/PEL 
approach r a is conveniently defined via r a = J dtSo(i) whereSo (t) denotes the probability 
that, starting form a randomly chosen time, the system has not performed a relaxation 
process [27J. It can be interpreted as the persistence time distribution [31] . Straightforward 
calculations yields r a = (r 2 ) /2(r) [221 E2]. So(t) is shown in Fig. [2] for all considered 
temperatures. For MD trajectories r a is typically extracted from the incoherent scattering 
function S(k.t). For a comparison of (ii), i.e. the MB hopping trajectory, with (iii) one 
has to choose a very large value of the wave vector k so that after one discrete hopping 
process full decorrelation is achieved. Choosing k = 400, one can see in Fig. [2] a very good 
agreement between Smb^ — 400, t) and S (t). For T = 0.6 it turns out that S (t) is slightly 
more non-exponential. One may speculate that this indicates the onset of anharmonic effects 
which at temperature T = 1 leads to a total breakdown of the PEL approach [331 134] . For 
the determination of r a we fit Smb with a stretched exponential function exp[— {t/To) KWW \- 
The a-relaxation time r a can then be calculated via r a = r (k c ) / (3Kww{k c )Y{l/ 0Kww(kc)) ■ 
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FIG. 2. Comparison of 5 (t) with S M B(k = 400, t) and S cont {k c = 15, t) for T = 0.45,0.5 and 0.6. 



As expected from the good agreement in Fig. [2] also the r a values of (ii) and (iii) agree 
well as shown in Fig. [T] More subtle is the quantitative comparison with (i) because the 
additional presence of vibrational and intra-MB processes give rise to additional decorrela- 
tion mechanisms for S con t(k,t). These aspects have to be worked carefully because for large 



system, see Sect. |III[ we are restricted to use the continuous MD trajectories. We start 
by fitting S con t(k,t) by a stretched exponential for times in the a-regime, yielding fitting 
parameters TQ 0nt (/c) and Pxwwi^)- Their fc-dependence for T = 0.5 is shown in Fig. [3j 
The same procedure is performed for the MB trajectory with the corresponding parameters 
r^ IB (k) and /S^ww^)- T o IB ( k ) and r^ oord (k) (calculated from MB respectively real space 
coordinates) correspond to each other up to a length scale of k ~ 10. For larger values of k 
one finds vibrational dynamics and T^ B (k) saturates, while TQ° ord (k) further decreases. The 
plateau values of Tq /b (/c) andf3^^r W (k) we are interested in can now be estimated by finding 
the value of k c where TQ OOTd {k c ) and Pkww^-c) matches with the MB plateau. We found 
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FIG. 3. k dependence of the fit parameter tq and @kww for the BMLJ65 at T = 0.5 for MB and 
real space coordinates. 

k c « 15. It is promising and non-trivial that for the same value of k both parameters can be 
recovered. k c is slightly temperature dependent, e.g. for T = 0.6 one gets k c m 14 ± 2. As 
shown in Fig. (2jS con t(k = k c ,t) agrees indeed very well with Smb^ = 400, t). Note that for 
this comparison the short time decay of S cont (k = k c ,t) has been scaled out. An equivalent 
scaling has been already used in Ref. [35] where continuous and IS trajectories had been 
compared. 

III. THE MACROSCOPIC SYSTEM 
A. General 



With increasing number of particles the properties of MB becomes less useful due to the 
following reasons: 1) For the minimum system consecutive MB transitions are uncorrelated, 
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because that during a transition the mobile particles are more or less equally distributed over 
the sample. In large systems, due to the dynamic heterogeneity, transitions are typically 
performed by the same mobile particles. Thus, the real space information about the location 
of a relaxation process is getting important. 2) Due to the many rearranging regions in the 
sample, the entire system performs a transition in any given fixed time interval, leading to a 
5-peaked waiting time distribution. 3) p e? (e) is narrowing when increasing the system size. 
As a consequence the relation between T and e is smeared out and the escape rate is much 
less correlated with the energy. 

As reported in [8], local waiting times remain a precise measure for the dynamics even 
in large systems. This means that CTRW results can also be applied in large systems. For 
the comparison with the a-relaxation time r a of the large system we have used S con t(k c ,t) 
since SMB(k,t) can no longer be defined (see discussion above). 

B. Evidence for coupling in a non-equilibrium configuration 

So far nothing is known about the size or shape of subsystems, furthermore the kind of 
coupling is unclear. But how can one identify coupling processes? When studying large 
systems one faces the problem of strong "dynamic noise" which makes it difficult to distin- 
guish between different coupling events. To minimize the dynamic noise, we prepared a very 
immobile configuration of a 520 particle system by copying a very stable N = 65 structure at 
T = 0.5. From this non-equilibrium structure we calculate the iso-configurational ensemble 
(IC) to study the distance dependence of possible coupling events between adjacent regions. 

At the beginning of each simulation two distinct behaviors can be identified: The a- 
particles are organized in a stable a-matrix which allows string like motion of b particles 
without changing the structure significantly. Hence we neglect b particles in the local event 
calculation. If an a-particle changes its position, it is either a subset of particles exchanging 
their positions inside the unaltered matrix or it reflects a significant local change of the 
matrix structure. The first process can be detected via local events but is irrelevant for the 
decay of the initial structure. The latter gives rise to define structural events: When an 
a-particle performs a local event at time t, we calculate the time averaged coordinates of 
the tagged particle and the first shell particles at t ± At with At m r a . To determine if the 
structure of the first shell changes significantly during the exchange process, we calculate 
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FIG. 4. Number of first events in the distance d normalized by the average number of particles 
one would expect from the radial distribution function. 

the squared displacement MSD = Xli( r *(^ + — r i'(^ — ^t)) 2 °f particle positions before 
and after the central event. The index il corresponds to the particle residing closest to the 
initial position of particle i (allowing permutation). In equilibrium one obtains (MSD) « 6, 
the distribution is very similar to a Gaussian with a variance close to 2. For the definition 
of a structural event we use a threshold of 2, above which we call the event structural. The 
precise choice of the threshold does not change the result. 

Before we present the MD results we first discuss possible effects in a small model system: 
Consider three independent systems with rates Tj = T. After the first (the left) system 
changes its state as a consequence of a relaxation process we now discuss the location of 
the next relaxation process. Here we concentrate on the middle and right system. More 
specifically we discuss the ratio (N m ) / (N r ) ((iVj): the number of next relaxation processes 
in system i) as well as w"}^^ / {ji^j (^ r i*2^ denoting the average waiting time between 

11 



the initial process of the left system and the next relaxation process in system i). Since 
both systems have the same rate one has (N m ) = (N r ) = 1 and {r^f) — V^F, yielding 
T \2 >Ss j I \ r i2/ = ^" Now we introduce a coupling mechanism where the initial relaxation 
enables the central system to acquire the rate T m = 2T whereas the rate T r of right system 
remains unchanged. Then one naturally has (N m ) / (N r ) = 2 and \Ti™' r ^ = 1/3T for 
both systems so that again \t[™^J / (^1^2^ = !• The situation changes if one assumes a 
distribution of rates p(T), e.g. T m = 2T and T m = 4T both with probability 0.5. In this 
case on has (N m ) I (N r ) = 3 and / (t$\ = 11/12 < 1. The ratio decreases with 

increasing width of p(T). 

For the BMLJ system this argument implies that as a consequence of a local coupling 
mechanism it is by far more likely that the second relaxation process occurs close to the initial 
one. Averaging over the IC ensemble would observe a significant distance ^-dependence 
N(d). The corresponding distance dependence of (t\^) gives information about the strength 
of the rate fluctuations. 

In Fig. [4] we show N(d) normalized by the average number of particles N av (d) one 
would expect from the radial distribution function g(r). The curve decays monotonically 
until a plateau value is reached for d > 2.5. As discussed for the model system this is a 
clear evidence for coupling processes. It would contradict the statistical case where no rate 
fluctuations are present. The curve roughly reaches a plateau value for d > 2.5 giving rise 
to the assumption that only particles within a sphere of radius r w 2.5, which is a little bit 
larger than the minimum system, are directly affected by the first relaxation event. Beyond 
this sphere, particles hardly recognize the first event, at first, and the changes of the rate 
become negligible. The discussion of 2 ) (d) can be found in Appendix [Cj 

The main contribution of the coupling is therefore between adjacent minimum systems 
so the coupling range can be restricted to such systems. These results reflect the underlying 
coupling processes and do not depend on the details of the definition of relaxation events. 



C. Physical picture 

To illustrate our physical picture of supercooled liquids we make use of the equilibrium 
distribution of rates p(T) introduced in (36] . Since thermodynamic properties [19] as well as 
the diffusivity [8] do not change upon increasing the system size (for iV > A" min ) we strongly 
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FIG. 5. Schematic sketch of the physical scenario: When the central system relaxes independently 
(active process, black arrow) the adjacent system may change its mobility without changing its 
state (passive process, dashed arrow). The role of q is discussed in the text. 



suggest that the local equilibrium distribution p(T) of the mobility does not change either. 
Therefore we assume that at any arbitrary time a macroscopic system can be decomposed 
into microscopic subregions of the size of the minimum system which can be described by 
the properties of their PEL, i.e. by values of e and T. When a subregion relaxes, this is 
called an active process. If one describes a real system by a decomposition into elementary 
systems, one has to include interactions between these subregions, which is what we call cou- 
pling. These interactions allow the subregions to change the rate while keeping their overall 
Boltzmann distribution p(T) (passive process), see Fig. |5} Now consider an immobile region 
adjoining to a rearranging mobile region: The adjacent relaxation enables rate fluctuations 
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in the slow sample, on average leading to a higher rate. This scenario directly explains the 
lack of slow regions in large model glass former compared to the minimum system jS]. The 
fluctuations themselves are interpreted as a coupling mechanism which leads to facilitation 
like dynamics, sometimes viewed as a hierarchical process [2]. The precise realization of 
passive processes will be discussed in the next sections. This kind of coarse graining and 
coupling is in principle captured in the kinetically constrained models, but two important 
differences remain: 1) All subsystems can always relax independently and 2) a subsystem 
already contains the complete macroscopic thermodynamics as also realized in the mosaic 
approach. 

D. The coupled landscape Model 

Since we have shown that the minimum system can be described by PEL parameters and 
the CTRW formalism and that localized coupling processes play a major role in supercooled 
liquids, we bring together both ingredients: We interpret a macroscopic glass former as a set 
of weakly coupled elementary systems which can be described by their PEL. The elementary 
systems (ES) are arranged on a square lattice, and their time evolution can be simulated 
as in [22]. After each active process all coupled adjacent systems are allowed to perform a 
passive process. We denote this approach as coupled landscape model (CLM). It is somehow 
a minimal model based on the PEL. It has to fulfill the condition that the thermodynamic 
properties as well as the diffusion coefficient D remain unchanged when increasing the system 
size. 

Dynamic coupling enables so called passive processes: Due to active processes of adjacent 
regions, the mobility (the rate) of an ES can also change without performing a relaxation 
process. The presence of passive processes must not change the equilibrium distribution 
p(e). Therefore the transition probability 7r(e ^ — > e new ) to move from state to state 
e new for a passive process has to fulfill 



The latter condition is needed for the probability interpretation. In this paper we focus on 
the most simple case n(e id — > e new ) = p(e new ), what we call Boltzmann coupling: With 
the probability q the new rate is simply chosen from the equilibrium distribution. Another 
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simple case is the Gaussian coupling where e new and e id are correlated, meaning that only 
small energy changes are allowed. With mean \i and variance a of the DOS one can define 
q(e id, ^new) = 1/V27rs 2 exp [— (e new — ae Q id — b) 2 /2s 2 ] with constants a = y/l — s 2 /a 2 and 
b = n(l — a). 

IV. RESULTS 

In this section we fix the coupling constant q and show that this model allows a non- 
trivial reproduction of the most important transport coefficients. This comparison is based 
on observables, which are well-defined in the MD simulation and the CLM. Here we estimate 
the coupling strength q based on the a-relaxation time r a . 

In previous work this kind of dynamical coupling was first estimated within a mean-field 
approach from the reduction of r a when going from N = 65 to iV = 130 particles by a 
linear expansion in q of the relaxation time r a (N) of the entire system [8]. Here, we follow 
a more general route which, first, is not based on linear expansion and, second, reflects the 
transition to macroscopic systems. For the simulations of the CLM we have used a 3 3 system 
with the parameters listed in table [T| 

A. a- Relaxation 

Again we will use finite size effects of r a to estimate q. We increase the system size 
from iVmin to N = 8320 and calculate T a (N) from the MD data. N = 8320 corresponds 
to the macroscopic limit for these temperatures. In the CLM one can model this scenario 
by adding additional ESs. If an ES corresponds to the smallest system in the MD, then 
one can directly compare the relative reduction of r a from an equilibrium simulation. For 
the dynamical interaction all systems sharing a boundary with the active can experience a 
passive process. 

The MD results show minor finite size effects for D [8]. In several papers this phenomenon 
is related to hydrodynamic interactions of the sample with its images due to the periodic 
boundary conditions [371 EE]- However, these hydrodynamic interactions are not captured 
by the model (D remains constant under coupling). To compare MD and model data, we 
first have to remove the hydrodynamic effects by scaling D N and r a by D N /D min . The 
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FIG. 6. Relative reduction of r a vs. coupling probability q for different temperatures. The solid 
dots correspond to MD results (r Q (8320)/r a (65) = 0.32,0.11,0.08 for T = 0.6,0.5,0.45), the gray 
region corresponds to the error intervals. 

correction of r a for the largest systems studied is about one order of magnitude smaller 
than the observed finite size effect itself, so these hydrodynamic effects are negligible for this 
study. 

We determine q by the condition that the reduction of r a of the BMLJ system and 
the CLM exactly match. In Fig. [6] we show the results for three different temperatures. 
From this analysis we get a temperature dependent coupling constant q(T). In the chosen 
temperature interval, q decreases by a factor of roughly 3. For small and intermediate values 
of q a the empirical formula T a (q) = r Q (0)/(l + Cy/q) fits the data very well. We mention 
by passing that the reduction does not seem to be analytical in the q — > limit, unlike the 
linear dependence in q we used earlier in a mean-field approach [8]. 

The estimation of r a provides data for the temperature dependence of q. To answer 
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FIG. 7. Coupling constant q vs. average waiting time for different temperatures. The black solid 
line is a power law as guidance to the eye. 

the question how the significance of the facilitation procedure evolves with temperature 
we compared q with typical timescales in the system. In Fig. [7] a sublinear scaling with 
1/ (r) can be found, meaning, that with decreasing temperature the number of successful 
passive processes decreases. At the same time the heterogeneity of the elementary systems 
increases making a passive process more effective. Both mechanisms compete and we can 
not determine the total effect on the significance of facilitation. A similar behavior was also 
found in granular systems [I], where the number of facilitation processes decreases when 
increasing the density up to the granular glass transition. 

Since r a depends on the first two moments of tp and the changes of (r) due to the 
coupling mechanism are negligible by construction, the reduction of r a in Fig. [6] mainly 
reflects the strong decrease of (t 2 ). In Fig. [8] we show the violation of the Stokes-Einstein 
relation to demonstrate the influence of the coupling processes on the dynamics at different 
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temperature. The upper black solid line is the analytical result without coupling. The 
dashed curves correspond to CLM data with constant coupling strength for different values 
of q and the red line to CLM data with a temperature dependent q with an exponent of 
£ = 0.5. The decoupling of D and r a is obviously strongly reduced by the coupling and £ 
does not depend on q in a first order approximation. For low temperatures it reaches the 
value for large q of approximately 0.175. The prefactor shows a strong q dependence. If one 

1 /2 

assumes a temperature dependent q, for example q ~ (r) 1 , see below for the motivation 
for this choice, it is also possible to reach intermediate values for £. In literature one finds 
values between 0.05 <£< 1/3, also values around 0.5 are reported for polymers [3"§HITj ). 

In Fig. [9] one can see that the CLM results, r a and also the non-exponentiality parameter 
Pkww-i match very well with the BMLJ8230 at k c . At long times the relaxation of the 
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FIG. 8. Effect of the passive processes on the Stokes-Einstein relation: Black lines are calcu- 
lated with fixed coupling constants, the red one corresponds to a temperature dependent coupling 
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constant, oc (r) ' , the green points correspond to the MD data. 
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FIG. 9. Overview over S(k c ,t), S con t(k c ,t) (both scaled for optimal overlap) and So(t) from the 
CLM for T = 0.6, 0.5, 0.45 and N = 65, 8320. 

CLM is a little bit faster than in the MD, at shorter times we find the inverse relation. 
However, the mismatch is small and can be rationalized by possible heterogeneity of the 
coupling: For the simulation we have used a fixed coupling constant. One can imagine, 
that in the microscopic system one has a distribution of coupling constants. The presence 
of smaller q values immediately leads to slower relaxation at long times, while the presence 
of larger coupling constants will increase the decay of the relaxation function at short and 
intermediate times. 



B. Discussion 

First we want to check how different choices for the geometry of possible passive processes 
influence our results. Let N e s denote the number of affected systems. In the last section 
an active process triggers all neighboring ES (sharing a boundary), i.e. iV eff = 2d. Other 
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realizations are possible, e.g. the elastic case: In [9] the authors studied finite size effects of 
the mechanical loss in thin films of a metallic glass and found that the loss vanishes below 
a certain thickness of the film. In their view the Eshelby stress field around a plastic zone 
is up to three times larger than the excitation itself. In dielectric relaxation experiments 
on glycerol it was found that at low temperatures the relaxation becomes nonlocal and r a 
decreases with system size [TT]. The latter effect was not observed in the CLM. However, 
these effects give rise to some long-range elastic interaction between excitations: An ES can 
be facilitated with the probability p(r) ~ r~ n for r < r max . Hence, every system contributes 
with r~ n to the number of effective interacting systems N e g. This value can be calculated 
as iVeff = Yli r ~ n - A s the coupling constant q only describes the probability to perform 
a passive process for one system, it is useful to introduce the effective coupling strength 



q e s = N e sq. Fig. [10] shows r a vs. g e ff for different interaction lengths r max and illustrates 
the importance of q e s as an effective parameter. All curves almost lie on a master curve, so 
that the distance dependence of the coupling effects only play a very minor role. 

One of the most important results is the temperature dependence of q. How can this 
be interpreted on the level of a single MB transition? One can think of the following 
physical scenario: Every state has an intrinsic resistance against being facilitated leading 
to an energy-dependent coupling constant q(e): the lower the energy, the more stable is 
the MB. This additional factor has, of course, to be considered when choosing the new 
state to obtain detailed balance. A general e dependence cannot be handled analytically, 
but for q(e) = qT(e) one can use the transition probability tt oc Tip which generates the 
correct statistics. If we now calculate q for different temperatures so that T a (g)/r Q (0) exactly 
matches with r a (q)/r a (0) we find some temperature dependence for q. Effects along this 
line may thus rationalize the observed temperature dependence. 



C. Conclusion 

In a first step we have analyzed in detail how the dynamics of a small BMLJ system with 
65 particles can be expressed in terms of PEL parameters. A key step in this endeavor is 
the fragmentation of the configuration space into MBs, each of which is characterized by 
an energy. Furthermore, in a simple PEL approach the rate to escape from a given MB is 
completely characterized by its energy and, more specifically, can be expressed by an effective 
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q e f f 

FIG. 10. Relative reduction of r a vs. coupling probability q e g for different neighbor geometries. 
Long range corresponds to r max = 5, intermediate range to r max = 2. 

barrier height and an entropic pref actor. Then it is possible to predict the dynamics for all 
temperatures. For this present case one additional aspect has to be taken into account: a 
given energy does not fully characterize the escape rate but the escape rate rather follows 
a relatively narrow distribution as expressed by a log-normal distribution. Intuitively, this 
reflects the fact that the BMLJ65 system can be interpreted as a superposition of approx. two 
elementary subsystems [26J. All PEL properties can solely be inferred from an appropriate 
analysis of the simulation data. 

As a consistency check one can estimate, e.g., the diffusivity or the structural relax- 
ation time. Indeed, an excellent agreement is found. Stated differently: the dynamics of 
the BMLJ65 system is very well understood in terms of the PEL. Of course, due to the 
discretization the PEL approach cannot resolve properties of the /3-relaxation. 

In principle one can perform the same MB discretization also for much larger systems. 
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However, in this case the total energy looses the tremendous information it has for small 
systems. Because a large system can be decomposed into (roughly) independent smaller 
systems the same total energy may result from many different realizations. As a consequence 
many different relaxation rates to escape from this MB are possible. Whereas for small 
systems the distribution of rates for a given energy is very small as compared to the total 
distribution of rates (as inferred from the different waiting times) this relation is inverted 
for large systems. Furthermore, for large system the mapping on the CTRW description is 
invalidated because of spatial relations. 

Therefore one has to complement the PEL approach of small systems by a new con- 
cept which takes into account the fact that different regions of the glass-forming system 
act somehow independently. Our CLM approach is guided by the observation of causal 
relations between different relaxation processes. Starting from an immobile non-equilibrium 
system these causal relations can indeed been identified. Via Monte Carlo simulations of 
the collection of BMLJ65 systems, interacting via the appropriately chosen coupling rules, 
relevant observables can be determined and compared with the properties of the large BMLJ 
systems. Interestingly, the precise choice of the coupling rule is not important for the prop- 
erties of D and r a . We carefully identified the wave vector for which the r Q -value has to 
be extracted in order to be compatible with the generic structural relaxation time from the 
CTRW approach. Comparison of the CLM with the actual system allows one to identify a 
coupling constant q which turns out to be temperature dependent. One may envisage that 
the resistance of a local region to change its relaxation rate as a consequence of close-by 
rearrangements becomes larger for lower local energies, then this T-dependence follows quite 
naturally. Most importantly, the variation of the whole shape of the incoherent scattering 
function S cont (k,t) when going from small to large systems can be reproduced by the CLM 
approach, i.e. by adjusting a single parameter. As a consequence the PEL parameters, 
defined for the small system, also determine the dynamics of large systems if supplemented 
by the coupling constant q. 

The CLM can be interpreted as a minimum approach to incorporate the dynamic cou- 
pling effects which is compatible with the key observation that the thermodynamics as well 
as the diffusivity only show very small finite size effect. On a qualitative level our model 
resembles the facilitation approach since immobile regions are typically rendered mobile by 
the properties of nearby regions and the elementary system size is temperature indepen- 
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dent. There is, however, one important difference. Whereas in the facilitation approach the 
elementary systems are just the spins and therefore do not contain any relevant thermody- 
namic or dynamic information, the elementary building blocks in the present case are small 
systems which already contain the (nearly) complete information about the thermodynamics 
and the diffusivity. 

In the literature one finds various arguments explaining the observed finite size effect 
that increasing the system size leads to a decrease of the relaxation time. However, in a 
comprehensive study Berthier et al. showed that this behavior is not universal [32]. If the 
relaxation time increases with increasing system size this has been related to an activated 
mechanism (defect diffusion). In contrast, the opposite behavior has been explained in terms 
of a mode-coupling like mechanisms where the cooperative relaxation occurs via unstable 
modes. Going to small systems unstable modes may disappear. This explanation is very 
different to the present PEL approach because here the relaxation mechanism is described 
by a trapping-type picture rather than by unstable modes. 

Karmakar et al. found a remarkable correlation of the relaxation time with the config- 
urational entropy [33], supporting the mosaic approach. The data suggest that a system 
containing approx. 1000 particles can serve as unit system in terms of the configurational 
entropy which is much larger than the building blocks in the CLM. However, it remains 
open why already for much smaller systems sizes the diffusivity displays macroscopic behav- 
ior. In [33] the authors explain the finite size scaling behavior of r a with the existence of a 
static length scale £(T) and an entropic argument: If the system size L is smaller than £(T) 
relaxation processes have to occur on the scale L with a smaller degeneracy factor. In larger 
systems the degeneracy factor will grow, leading to a lower free energy barrier. In spirit 
this argument resembles the mode-coupling mechanism from [32]. For the Kob- Andersen 
system £(T) changes by roughly 30% within the chosen temperature interval and our data 
suggest that a minimum system with around 65 particles is large enough to reproduce the 
thermodynamics as well as the diffusivity in this temperature interval. Since the minimum 
system does not show a significant temperature dependence, we can only speculate that 
the relevant static length scales for all temperatures are already captured by the minimum 
system. 

This entropy-related argument, and the CLM differ in one important point: While the 
earlier one is based on the free energy barrier height at a fixed time, the finite size effect in the 
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CLM is a consequence of fluctuating barriers (and corresponding jump rates) over a certain 
time interval induced by the facilitation. Since the free barrier from [H] also determines 
the self diffusion it remains unclear how the different scaling behavior of D and r a can be 
reconciled with this approach. In the CLM, at a fixed time the distribution of rates agrees 
with the equilibrium distribution of the minimum system. From this distribution directly 
follows the lack of finite size effects for the thermodynamics and the diffusivity, while the 
fluctuations effect r a exclusively. 

It may be interesting to study in future work whether it is also possible to specifically 
describe the viscosity of large systems via a coupling of small systems. The observed strong 
correlations between structural relaxation time and viscosity suggest that a similar approach 
should be possible. 
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Appendix A: Analytical results 

Even in the generalized case it is possible to calculate (T n (e)) = J dep eq (e)T n (e) an- 
alytically as needed for the calculation of D and r a . Solving the Gaussian integrals one 
obtains 
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the transport coefficients then read 
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FIG. 11. The different waiting times (for the definition see text) for a MSD threshold of 2 vs. the 
corresponding distance. 

Appendix B: How to handle a system with A < 1 

When modeling a BMLJ65 with e.g. two independent PEL systems one faces a problem: 
In the MD one can only access r 2 (e). If one calculates the entire rate T^(e) of two PEL 
systems with Ti(e) for the model system from [22], the superposition leads to an additional 
temperature dependent factor F(T) (from the thermodynamic distribution of energies). In 
the general case it is even possible to end up in an additional energy dependence, making it 
impossible to extract Ti(e) from the MD data. 

When one uses a single system with A < 1 one hast to take into account that the 
persistence time distribution is no longer exponential and has to be determined numerically 
when simulating the system. Furthermore, (T 2 (e)) / (T(e)) 2 > 2 as obtained in the MD 
simulation has to be included separately by the distribution e) (see text). Another 
result is that the invariance of (r) under passive processes is no longer valid since (r(e)) = 
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1/ (r(e)) does not hold for <7r > 0. However, this effect is very small and even for high 
coupling strength smaller than w 4% (at T = 0.45). 



Appendix C: Waiting times of the non-equilibrium ICE 

As discussed in the text a possibility to gain information about the rate fluctuations of the 
coupling is the waiting time (t 12 ) (d) between the first two structural events in dependence 
on their distance d. In the case that subsequent relaxation processes are uncorrelated or the 
coupling mechanism uses a fixed rate one would expect (ti j2 ) (d) = (to,i). In the presence 
of coupling with a rate distribution, small d will display small (ti^) (d). The results are 



shown in Fig. 11 We found (ti^) increases monotonically with growing distance d up to 
half of the cell length. Very small values of (ti 2) are seen for d < 2.5, followed by a large 
jump in (t 12 ). These findings again contradict the statistical case, but since we are in a 
non-equilibrium configuration we cannot extract information about the strength of the rate 
fluctuations. The waiting time 7"o,i for the first structural event is also included in the figure. 
The matching of both observables at large d is to our knowledge highly nontrivial. 
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